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ABSTRACT 

We present numerical simulations on propagation of Ultra-High Energy Cosmic Rays (UHECRs) above 10 19 
eV in a structured extragalactic magnetic field (EGMF) and simulate their arrival distributions at the earth. We 
use the IRAS PSCz catalogue in order to construct a model of the EGMF and source models of UHECRs, both of 
which reproduce the local structures observed around the Milky Way. We also consider modifications of UHECR 
arrival directions by the galactic magnetic field. We follow an inverse process of their propagation from the 
earth and record the trajectories. This enables us to calculate only trajectories of UHECRs arriving at the earth, 
which saves the CPU time. From these trajectories and our source models, we construct arrival distributions of 
UHECRs and calculate the harmonic amplitudes and the two point correlation functions of them. We estimate 
number density of sources which reproduces the Akeno Ground Air Shower Array (AGAS A) observation best. As 
a result, we find that the most appropriate number density of the sources is ~ 5 x 10~ 6 Mpc" 3 . This constrains the 
source candidates of UHECRs. We also demonstrate skymaps of their arrival distribution with the event number 
expected by future experiments and examine how the EGMF affects their arrival distribution. A main result is 
diffusion of clustering events which are obtained from calculations in the absence of the EGMF. This tendency 
allows us to reproduce the observed two point correlation function better. 

Subject headings: cosmic rays — methods: numerical — IGM: magnetic fields — galaxies: general — 
large-scale structure of the universe 



1. INTRODUCTION 

The nature of Ultra-High Energy Cosmic Rays (UHECRs), 
which are particles of energy above 10 19 eV, is poorly known. 
This is one of the most challenging problems of modern astro- 
physics. 

One of problems about UHECRs is what their origin is. The 
two scenarios of their origin are suggested, which are called 
bottom-up and top-down ones. On the one hand, bottom-up 
scenarios assume some astrophysical phenomena as their ori- 
gin. UHECRs are thought to be of extragalactic origin since 
the gyroradii of UHECRs above 10 19 eV exceed thickness of 
, our galaxy. From this fact and the Hillas plot (Hillas 1984), 
probable candidates of UHECR origins are active galactic nu- 
clei (AGNs), gamma-ray bursts (GRBs) and colliding galaxies. 
Theoretically, this scenario predicts the GZK cutoff of the en- 
ergy spectrum of UHECRs (Greisen 1966; Zatsepin & Kuz'min 
1966) since their source candidates are located at far distances. 
UHE protons with energy above ~4x 10 19 eV interact with 
the cosmic microwave background (CMB) and lose large frac- 
tion (~ 20%) of their energy per interaction by photopion pro- 
duction (Berezinsky & Grigorieva 1988; Yoshida & Teshima 
1993). The mean free path of UHE protons through the CMB 
field is ~ 10 Mpc at 10 20 eV. Thus the energy spectrum at the 
Earth should have a cutoff around £~8x 10 19 eV. This spec- 
tral cutoff is called the GZK cutoff. But there is a observa- 
tional disagreement of the energy spectra between the Akeno 
Giant Air Shower Array (AGAS A), which does not observe the 
GZK cutoff (Takeda et al. 1998), and the High Resolution Fly's 
Eye (HiRes; Wilkinson et al. 1999), which does it (Abu-Zayyad 
et al. 2004). This discrepancy between the two experiments 
remains being one of open questions in astroparticle physics. 
On the other hand, top-down scenarios assume some processes 



based on new physics beyond the standard model of the particle 
physics (see a review Bhattacharjee & Sigl (2000)). 

Another problem is arrival distribution of UHECRs. The 
AGASA reported that there is no statistically significant large 
scale anisotropy in the observed arrival distribution of UHE- 
CRs above 10 19 eV (Takeda et al. 1999). This fact points 
out that sources of UHECRs are distributed isotropically, but 
isotropic distribution of sources cannot reproduce the small- 
scale anisotropy reported by the AGASA (Takeda et al. 1999, 
2001). A model of UHECR origin are constrained by their abil- 
ity to reproduce such observed arrival distribution of UHECRs. 
On the other hand, the HiRes experiment indicates that there 
is no statistically significant small-scale anisotropy (Abbasi et 
al. 2004; Farrar et al. 2004). However Yoshiguchi et al. (2004) 
concluded this discrepancy between the two observations is not 
statistically significant at present. This problem is left for fu- 
ture investigation by new experiments such as the Pierre Auger 
Observatory. 

To obtain information on UHECR origin, we need to calcu- 
late their arrival distribution using some kinds of source mod- 
els. To do so, we have to simulate propagation of UHECRs in 
the intergalactic space, where the Extragalactic Magnetic Field 
(EGMF) plays important roles since we assume that UHECRs 
are protons in this paper. Yoshiguchi et al. (2003a) calculated 
propagation of UHE protons in an uniform turbulence of mag- 
netic field with the Kolmogorov spectrum. But such magnetic 
field is not realistic since the EGMF is expected to reflect the 
large scale structure of the universe. 

In recent years, several groups started to develop physically 
more realistic models of the EGMF based on numerical simula- 
tions of large scale structure formation. Sigl, Miniati, & Ensslin 
(2003, 2004) used a structured EGMF model which is gener- 
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FIG. 1 . — The IRAS PSCz galaxies within 100 Mpc on the galactic coordinate. Strong infrared luminosity on the galactic plane interrupts observation of galaxies 
near the galactic plane. 



ated by their large scale structure simulations, taking magnetic 
fields into account. But their model does not reproduce the local 
structures actually observed around the Milky Way. This causes 
the ambiguity in the choice of observer position. In addition 
to this, the calculated arrival distribution of UHECRs does not 
correspond to the one expected at the earth. 

An important step on modeling the magnetic structure of the 
local universe is performed by Dolag et al. (2005). They con- 
strain the initial conditions for the dark matter density fluctua- 
tions to reproduce the local structures. This allows us to remove 
the ambiguity in the choice of observer position, and to ob- 
tain the simulated skymaps of expected UHE proton deflections 
in the magnetic large-scale structure around our galaxy. How- 
ever, they did not calculate the arrival distribution of UHECRs, 
and thus could not obtain the information on source distribu- 
tion which reproduces the AGASA observation. And also, the 
effects of the Galactic Magnetic Field (GMF) are not consid- 
ered in both Sigl, Miniati, & Ensslin (2003, 2004) and Dolag et 
al. (2005). But recently it has also been shown that the GMF 
affects the arrival distribution of UHECRs (Alvarez-Muniz, En- 
gel & Stanev 2002; Yoshiguchi et al. 2003b). Thus we cannot 
neglect modifications of arrival directions of UHECRs by the 
GMF when we simulate their arrival distribution. 

In this work, we calculate propagation of UHECRs taking 
both the EGMF and the GMF into account, and simulate the 
arrival distribution of UHECRs. We constrain source num- 
ber density of UHECRs by comparison of the results with the 
AGASA observation. We generate the magnetic structure of 
the local universe by our original method (section 3.1) from the 
IRAS PSCz catalogue of galaxies (Saunders et al. 2000). We 
also construct our source models of UHECRs from this cata- 
logue. As our GMF model, we adopt a bisymmetric spiral field 
(BSS) model (Alvarez-Muniz, Engel & Stanev 2002) just like 
Yoshiguchi et al. (2003b). 

In order to simulate the arrival distribution of UHECRs, we 



apply a method developed in previous works. (Fliickiger et 
al. 1991; Bieber, Evenson, & Lin 1992; Stanev 1997; Medina- 
Tanco 1999; Yoshiguchi et al. 2003b). We numerically calcu- 
late an inverse process of propagation of UHE protons, which 
reach the earth, and record their trajectories in our Galaxy and 
the intergalactic space. In other words, we inject UHECRs 
from the earth isotropically whose charges are taken as -1. We 
then select some of them according to a given source distri- 
bution. (Detailed explanation is given in the section 4.3) The 
expected arrival distribution can be obtained by mapping the 
velocity directions of the selected trajectories at the earth. The 
validity of this method is supported by the Liouville's theorem. 
This method enables us to save the CPU time effectively since 
we calculate only trajectories of UHE protons which reach the 
earth. A method for this process is explained in section 4. 

The outline of this paper is as follows. In section 2, we ex- 
plain the IRAS PSCz Catalogue and construct our sample of 
galaxies. In section 3, we introduce our model of the EGMF 
and the GMF. We explain our numerical methods for calculat- 
ing arrival distribution of UHECRs and statistical quantities in 
section 4. Then in section 5, we estimate the most appropriate 
number density of source of UHECRs and compare the statisti- 
cal quantities calculated from this source model with the EGMF 
to those calculated without the EGMF. We also demonstrate 
skymaps of the arrival distribution of UHECRs. In section 6, 
we summarize the main results. 

2. A SAMPLE OF GALAXIES 

In order to calculate propagation of UHECRs considering the 
local structures actually observed around the Milky Way, we 
use the IRAS PSCz Catalogue (Saunders et al. 2000) for con- 
struction of our EGMF model and UHECR source models. 

We used the ORS sample of galaxies to construct our 
UHECR source models in our previous work (Yoshiguchi et 
al. 2003a). The ORS sample has better completeness on nearby 
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galaxies than the IRAS PSCz catalogue. Thus we used the ORS 
sample to construct our source model since we were interested 
in the local sources in order to research nature of UHECRs 
above 4 x 10 19 eV. As a result, we obtained that UHECR source 
number density is ~ 10~ 6 Mpc~ 3 . Then we investigated that this 
source model also explained arrival distribution of UHECRs 
above 10 19 eV though the EGMF was neglected (Yoshiguchi 
et al. 2003b). In this work, we use galaxy sample to construct 
not only a source model of UHECRs but also a model of the 
EGMF, which reflect the large scale structure of the universe. 
This requests large sky coverage of the galaxy sample. Thus 
we adopt not the ORS galaxy catalogue but the IRAS PSCz 
Catalogue. 

The IRAS PSCz catalogue consists of 14677 galaxies with 
redshift and infrared fluxes > 0.6 Jy, and covers about 84% 
of the sky. We assume de Sitter universe with Q m = 0.3, £1 a = 
0.7, H = 71 km s" 1 Mpc" 1 in order to calculate distance of each 
galaxy. We show distribution of the IRAS PSCz galaxies in 
Figure 1. 

However, there are two problems when we use the IRAS 
PSCz Catalogue. One is that it is impossible to observe dark 
galaxies at far-off distance, (the selection effect, see figure 2). 
The other is that this catalogue has the zone of avoidance (the 
mask), where the IRAS PSCz Survey did not observe galax- 
ies. Our previous works (Yoshiguchi et al. 2003a,b), which 
uses the ORS sample (Santiago et al. 1995), also have these 
problems. Using luminosity function of the IRAS galaxies, we 
correct these absence of galaxies in the same manner with our 
previous studies. 

We use the luminosity function of the IRAS PSCz galaxies 
(Takeuchi et al. 2003), 



*CL) = exp \-±- 2 



log 



i + A 
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Here L„ = (4.34 ±0.86) x 1O 8 /T 2 [L ], a = 1.23 ±0.04, a = 
0.724 ±0.010,0* = (2.34 ± 0.30) x 10- 2 /z 3 Mpc 3 . Using this lu- 
minosity function, we define the selection function as, 



cj>{r) : 
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where L m i n (r) is minimum luminosity of galaxies which are ob- 
servable at a distance r. Therefore, 4>{r) represents fraction of 
all galaxies that are observable at each distance. 

First of all, we correct the selection effect. For each of the 
IRAS galaxies, we add galaxies that are not included in the 
IRAS sample. The number of added galaxies can be obtained 
by using the selection function. The positions of added galax- 
ies are determined according to the Gaussian distribution whose 
mean is the position of the original IRAS galaxy and whose root 
mean square is l(r). Here we define l(r) as a mean distance be- 
tween the original IRAS galaxies at distance r, thus 
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where n(r) is number density of the original IRAS galaxies. 
Then, 
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The luminosities of added galaxies are randomly assigned 
so that their distribution of luminosity is consistent with the 



luminosity function. This method can complement the IRAS 
galaxies without spoiling the observed structure of galaxy dis- 
tribution. 




FIG. 2. — Infrared luminosity (60/im) of the IRAS PSCz galaxies plotted as 
a function of distance from the Earth. Unit of luminosity is L© . 

Next, we add galaxies in the mask. We assume that galaxy 
distribution in this region is homogeneous and number density 
of galaxies is n{r)4>{r)~ l . These luminosities are random but 
their distribution is consistent with the luminosity function. Our 
galaxy sample after these corrections is shown in Figure 3. 

In this work, we use only the IRAS galaxies within 100 Mpc. 
We assume that source distribution at r > 100 Mpc is isotropic 
and uniform, and that their number density is equal to that 
within 100 Mpc. We neglect cosmological evolution of num- 
ber density of galaxies. Thus our sample of galaxies reflects 
the observed local structure within 100 Mpc. We use this sam- 
ple of galaxies to construct a model of EGMF and source model 
of UHECRs. 



3. A MODEL OF MAGNETIC FIELD 
3.1. Extragalactic Magnetic Field 

The EGMF are little known theoretically and observation- 
ally. Theoretically, several large scale structure simulations 
with magnetic field have been performed (Dolag et al. 2005; 
Sigl, Miniati, & Ensslin 2004). Roughly speaking, their results 
are that the strength of magnetic field traces baryon density. A 
model of Sigl, Miniati, & Ensslin (2004) do not reflect the local 
structures actually observed around the Milky Way, while one 
of Dolag et al. (2005) reflects these local structures. To com- 
pare model predictions of UHECR arrival distribution with the 
observed one, it is important to generate the magnetic structure 
around our galaxy. Therefore, we present a model of the EGMF 
reflecting the local structures of the universe well. 

The EGMF mainly exists in clusters of galaxy or around 
galaxies. From this standpoint, we present a realistic model 
of the EGMF. We assume that magnetic field results from the 
amplification of weak seed fields. In a simulation of evolution 
of cluster of galaxy (Dolag et al. 2002), the average magnetic 
field strength in the cluster is amplified as expected from com- 
pression alone (|B| oc p 2//3 , where p is density of matter). We 
adopt this conclusion 



ocp 2 / 3 ocp L 2 / 3 , 



(5) 
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FIG. 3. — Distribution of galaxies after our correction for the selection effect and existance of the mask. This distribution do not spoil the observed structure of 
galaxy distribution. 



where pi is luminosity density explained in the next paragraph. 
In equation 5, we assume that the luminosity density on each 
position is proportional to density of the gas. With these as- 
sumption, we construct a model of the EGMF as follows. 

First, we cover the universe with cubes of side l c , which is 
correlation length of the EGMF. We adopt l c = IMpc from our 
previous work (Yoshiguchi et al. 2003a). Second, we sum lu- 
minosities of galaxies in our sample which exist in each cube. 
We call this summed value luminosity density. Finally, It is 
assumed that the magnetic field in each cube is represented as 
the Gaussian random field with zero mean and has a power-law 
spectrum 

(B(t)B*(t)j oc k" H for 2tt/I c <k< 2ir/l cutl 

B(T)B*(T)^> = otherwise, (6) 

where Z cut is a numerical cutoff scale. Physically, one expects 
'cut *C k, but we set l cut = 1 /8 x l c in order to save the CPU time. 
We use n H = -H/3, corresponding to the Kolmogorov spec- 
trum since the Faraday rotation map reveals that the clusters' 
magnetic fields are turbulent with the Kolmogorov spectrum 
over at least one order of magnitude of the wavenumber(Vogt 
& Ensslin 2004). 

Next we consider a normalization of the EGMF. Most ob- 
servations suggest that clusters of galaxy have magnetic field 
whose strength is from 0.1 /xG to a few /iG (see a review Vallee 
(2004)). On the one hand, the Faraday rotation measurements 
of polarized radio sources placed within cluster of galaxies pro- 
vide some evidence for the presence of stronger intracluster 
magnetic field (ICMF), in the range of a few /iG (Taylor et al. 
2001; Vogt et al. 2003). On the other hand, observations of 
hard X-ray emission from cluster of galaxy implies that an av- 
erage ICMF strength within the emitting volume is 0.2-0.4 /iG 
(Fusco-Femiano et al. 1999; Rephaeli et al. 1999). 
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FIG. 4. — Magnetic field strength along three fiducial lines through the Virgo 
cluster, the Perseus cluster and the Coma cluster within 100 Mpc. 

In this work, we now set a normalization of its strength 
~ 0.4/iG in a cube where is the center of the Virgo cluster. 
In order to compare our model with Dolag et al. (2005), we 
show magnetic field strength within 100 Mpc along three fidu- 
cial lines through the Virgo cluster, the Perseus cluster and the 
Coma cluster in figure 4. In figure 5, we also show deflection 
maps when protons propagate through our EGMF from the dis- 
tance of 100 Mpc. The deflection angle by uniform turbulent 
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magnetic field is given by 
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Here r is distance of propagation. In the way similar to Dolag 
et al. (2005), we assume that proton trajectory makes a random 
walk through each cell since our EGMF in each cube of side l c 
is the turbulence. These maps can be compared to figure 13 and 
14 in Dolag et al. (2005). 

We use our sample of galaxies to construct a model of the 
EGMF only within 100 Mpc. At a distance above 100 Mpc, 
we treat the EGMF as an uniform turbulence of magnetic field 
with the same spectrum and \B\ = 1 nG since the IRAS PSCz 
catalogue is poorly covered at r > 100 Mpc. 

3.2. Galactic Magnetic Field 

In this study, we adopt the GMF model used by Alvarez- 
Muniz, Engel & Stanev (2002), which is composed of the spiral 
and the dipole field. We briefly explain this GMF model below. 

Faraday rotation measurements indicate that the GMF in the 
disk of the Galaxy has a spiral structure with field reversals at 
the optical Galactic arms (Beck 2001). We use a bisymmetric 
spiral field (BSS) model, which is favored from recent work 
(Han, Manchester, & Qiao 1999; Han 2001). The Solar Sys- 
tem is located at a distance rn = R® = 8.5 kpc from the center 
of the Galaxy in the Galactic plane. The local regular mag- 
netic field in the vicinity of the Solar System is assumed to be 
#Soiar ~ 1-5 P-G m tne direction / = 90° +p where the pitch angle 
is p = -10° (Han & Qiao 1994). 

In the polar coordinates (n | , <j>), the strength of the spiral field 
in the Galactic plane is given by 



B(.r\\,<t>) = Bo 



r \\ 



cos 



-/31n 



(8) 



where B Q = 4.4 fiG, r = 10.55 kpc and (3 = 1/tanp = -5.67. 
The field decreases with Galactocentric distance as 1/r^ and it 
becomes zero for ni > 20 kpc. In the region around the Galac- 
tic center (n < 4 kpc) the field is highly uncertain, and thus 
assumed to be constant and equal to its value at rj| =4 kpc. 

The spiral field strengths above and below the Galactic plane 
are taken to decrease exponentially with two scale heights 
(Stanev 1997), 



| B(ni ,<A,z)| = |B(ni^)li exp(f)exp(=M) 



exp(-|z|) : 



\z\ < 0.5 kpc 
\z\ > 0.5 kpc 
(9) 

where the factor exp(-3/8) makes the field continuous on z. 
The BSS spiral field we use is of even parity, that is, the field 
direction is preserved at disk crossing. 

Observations show that the field in the Galactic halo is much 
weaker than that in the disk. In this work we assume that the 
regular field corresponds to a AO dipole field as suggested in 
(Han 2002). In spherical coordinates (r,9,(p), the (x,y,z) com- 
ponents of the halo field are given by: 

B x = -3 sin 9 cos 9 cos (p/r 3 

By = -3 sin#cos#sin(^/r 3 (10) 

B 7 = nc (1- 3 cos 2 <9)/r 3 
where — 184.2 /iGkpc 3 is the magnetic moment of the 
Galactic dipole. The dipole field is very strong in the central 



region of the Galaxy, but is only 0.3 /iG in the vicinity of the 
Solar system, directed toward the North Galactic Pole. 

There is a significant turbulent component, B ran dom, of the 
Galactic magnetic field. Its field strength is difficult to mea- 
sure but results found in literature are in the range of B ran dom = 
0.5 .. . 2B mg (Beck 2001). However, we neglect the random field 
through the paper. Possible dependence of the results on the 
random field is discussed in the section 5. 



4. NUMERICAL METHOD 
4.1. Method of Calculation for Propagation of UHECRs 

We numerically calculate an inverse process of propagation 
of UHE protons arriving at the earth in the intergalactic space. 
This method explained below is an expansion of many previous 
works on their propagation in the Galactic space (Fliickiger et 
al. 1991; Bieber, Evenson, & Lin 1992; Stanev 1997; Medina- 
Tanco 1999; Yoshiguchi et al. 2003b). 

We already performed numerical simulations for UHECR 
propagation in the GMF in Yoshiguchi et al. (2003b). In the 
paper, we injected UHECRs from the earth isotropically, and 
recorded these trajectories until they reached a sphere of radius 
40 kpc centered at the Galactic center. The charge of UHE- 
CRs was taken as -1 because we followed propagation of UHE 
protons backward. These UHECRs were injected with spectral 
index of -2.7, which was similar to the observed one. Note that 
this is not the energy spectrum injected at extragalactic sources. 

In this study, we expand these trajectories to the extragalactic 
space. In other words, there are our initial positions of UHE- 
CRs on a sphere of radius 40 kpc centered at the Galactic center, 
which are the result of our previous work. The trajectories are 
followed until their distance from the Galaxy reaches 1 Gpc or 
their time for propagation reaches the age of the universe or 
their energies reach 10 25 eV. Of course, we set the charge of 
UHECRs to be -1. 

In the extragalactic space, we have to consider not only the 
deflections due to the EGMF but also the energy loss processes 
(Berezinsky & Grigorieva 1988; Yoshida & Teshima 1993). 
UHE protons below 4 x 10 19 eV lose their energies mainly due 
to adiabatic energy losses and pair production in collision with 
the cosmic microwave background (CMB). At the higher ener- 
gies the photopion production with the CMB becomes essen- 
tial (Detail explanation is given below). Though we assume 
that UHECRs are protons in this work, we should also add the 
photo-disintegration if we assume UHECRs to be nuclei. We 
treat all these energy loss processes as continuous processes. 
Note that energies of UHECRs increase during propagation be- 
cause we follow their inverse processes. 

The adiabatic energy loss is the effect of the expanding uni- 
verse. This energy loss is written as 



^=--E = -H Q [n m (i + zf+n A ] E. 

dt a 1 J 



(11) 



As mentioned in the section 2, the cosmological parameters 
used in this calculation are fi m = 0.3, Oa = 0.7, and Hq = 
71 kms" 1 Mpc" 1 . 

The pair production due to collisions with the CMB can be 
treated as a continuous process which has small inelasticity 
(<~ 10~ 3 ). We adopt the analytical fit functions given by Chodor- 
owski, Zdziarske, & Sikora (1992) to calculate the energy loss 
rate on isotropic photons. 
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UHE protons above ~4x 10 19 eV lose a large fraction of 
their energy (~ 20% at every reaction) in the photopion pro- 
duction with the CMB. We treated this process as a stochastic 
process in previous work. But in this study, we cannot treat this 
as a stochastic process because we calculate the inverse pro- 
cess. Berezinsky, Gazizov & Grigorieva (2002) and Aloisio & 
Berezinsky (2004) showed that the energy spectrum which is 
calculated with a continuous process of the photopion produc- 
tion is consistent with one calculated with a stochastic process 
if a distance between the earth and sources of UHECR is more 
than 30 Mpc. Thus, we can adopt a continuous energy loss pro- 
cess since our source model (explained below) almost satisfy 
that condition about source distance. We use the energy loss 
length which is calculated by simulating the photopion produc- 
tion with the event generater SOPHIA (Mucke et al. 2000). 

4.2. Source Distribution 

We construct source models of UHECRs from our sample 
of galaxies explained in the section 3. The number density of 
UHECR sources is taken as our model parameter. For a given 
number density, we randomly select galaxies from our sam- 
ple with probability proportional to absolute luminosity of each 
galaxy. We then estimate the source number density which re- 
produces the observed arrival distribution of UHECRs, by cal- 
culating the harmonic amplitude and the two point correlation 
function of arrival distribution of UHECRs as a function of the 
source number density. 

4.3. Calculation of the UHECR Arrival Distribution 

In this subsection, we explain the method of construction of 
UHECR arrival distribution at the earth. 

We calculate 500,000 trajectories of UHE protons in the 
EGMF, using our method explained in Section 4.1, and record 
them. With our source models, we calculate a factor for each 
trajectory, which represents a relative probability that j th pro- 
ton reaches the Earth, 



■Pselec (£,./') OC 



dN/dE g (d Lj ,Ei}dE g 



iX+Zij)di/ 



-2.7 



dE 



(12) 



Here, i labels sources on each trajectory. Zij, dij and Ly is 
a redshift, a distance and luminosity of each source which is 
passed by j th proton respectively. dN /dE(dij,Ei) is the energy 
spectrum of protons at a source of distance d, j. £, is energy of 
proton at i th source. E g = E g (E,d) is the energy of cosmic 
ray at a source, which has the energy E at the earth. dE g /dE 



represents variation of shape of the energy spectrum through 
propagation. 

dE g /dE can be calculated in the case of rectilinear prop- 
agation (Berezinsky & Grigorieva 1988). But calculation of 
this is difficult in this study since protons injected from the 
sources which are located at the same distance have different 
path length due to the EGMF. That is, the only E g cannot be 
decided when E and d are given. We calculate this factor using 
our 500,000 trajectories of protons. 

Figure 6 shows a variation of shape of a monoenegetic spec- 
trum (E = 10 19 6 eV) at the earth for an example. The solid 
histogram is the spectrum at the earth. The dashed histogram 
and the dotted histogram are the spectra of UHECRs injected 
from the earth at 300 Mpc from us and 500 Mpc respectively. It 
is difficult to determine dE g from the figure since the spectra at 
far distances from the earth have large variances due to differ- 
ence of their path lengthes. In this case, we calculate dE g /dE 
as 



dE„ 
dE 



dN/dE(E) 



dN/dE g (E g *(E,d),d)' 



(13) 



where dN/dE(E) is the spectrum at the earth (oc E~ 2J ), 
dN /dE g (E g (E ,d),d) is a spectrum of UHECRs injected from 
the earth at a distance d and E g *(E,d) is averaged E g when E 
and d are given. 

We randomly select several trajectories according to these 
relative probabilities, so that the number of the selected trajec- 
tories is equal to the observed event number. The mapping of 
the velocity directions of each UHECR at the earth becomes the 
arrival distribution of UHECRs. The validity of this method is 
supported by the Liouville's theorem. 

If we have to select the same trajectory more than once in 
order to adjust the number of the selected trajectories, we gen- 
erate a new event whose arrival angle is calculated by adding 
a normally distributed deviate with zero mean and variance 
equal to the experimental resolution 2.8° (1.8°) for E > 10 19 
eV (4 x 10 19 eV) to the original arrival angle. 

We assume that UHECRs are protons injected with a power 
law spectrum in the range of 10 19 - 10 22 eV. We set this power 
law index 2.6 in order to fit the calculated energy spectrum 
to the one observed by the AGASA (Marco, Blasi, & Olinto 
2003). In other words, dN /dE(djj,E) oc E~ 2 6 where E is the 
energy of UHECR at the source. 
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Initial spectrum 
SpGCtrum at 300 Mpc 
SpGCtrum at 400 Mpc 



FIG. 6. — A variation of shape of a monoenegetic spectrum (E = 10 eV) 
at the earth through UHECR propagation. The solid histogram is the spectrum 
at the earth. The dashed histogram and the dotted histogram are the spectrum 
of UHECRs injected from the earth at 300 Mpc and at 500 Mpc from us re- 
spectively. The spectra at far distances from the earth have large variances due 
to difference of their path lengthes. 



4.4. Statistical Methods 

In this subsection, we explain the two statistical quantities, 
the harmonic analysis for large scale anisotropy (Hayashida et 
al. 1999) and the two point correlation function for small scale 
anisotropy. 

The harmonic analysis to the right ascension distribution of 
events is the conventional method to search for large scale 
anisotropy of cosmic ray arrival distribution. For a ground- 
based detector like the AGAS A, the almost uniform observation 
in right ascension is expected. The m-th harmonic amplitude r 
is determined by fitting the distribution of cosmic rays to a sine 
wave with period 2ir/m. For a sample of n measurements of 
phase, 0i, 02, ••;<f> n (P< 4>i < 2ir), it is expressed as 



where, a = 



■S^j cosm« 



r = {a 2 + b 2 ) 1 ' 2 
U, b= 2 -Y l " 



sin mo 



(14) 

We calculate the 



harmonic amplitude for m = 1,2 from a set of events generated 
by the method explained in the section 4.3. 

If events with total number n are uniformly distributed in 
right ascension, the chance probability of observing the ampli- 
tude > r is given by, 



P = exp(-fc), 



(15) 



where 



k = nr 2 /4. (16) 

The current AGASA 775 events above 10 19 eV is consistent 
with isotropic source distribution within 90 % confidence level 
(Takeda et al. 2001). We therefore compare the harmonic am- 
plitude for P = 0. 1 with the model prediction. 

The two point correlation function N(9) contains informa- 
tion on the small scale anisotropy. We start from a set of events 
generated from our simulation. For each event, we divide the 
sphere into concentric bins of angular size A9, and count the 
number of events falling into each bin. We then divide it by the 
solid angle of the corresponding bin, that is, 



N(0)-. 



1 



2tt|cos0-cos(0 + A0)| 



1 [ sr_1 ]< 



where <j> denotes the separation angle of the two events. A9 
is taken to be 1° in this analysis. The AGASA data shows 
correlation at small angle (~ 3°) with 2.3 (4.6) a significance 
of deviation from an isotropic distribution for E > 10 19 eV 
(£>4x 10 19 eV) (Takeda et al. 2001). 



5. RESULTS 

In this section, we present results of our simulations. In 
section 5.1, we constrain number density of UHECR sources 
from the observational results of the AGASA. Using our source 
model with this number density, we see how the EGMF af- 
fects the arrival distribution of UHECRs in section 5.2 and sec- 
tion 5.3. 



5.1. A constraint on source model of UHECRs 

In this subsection, we constrain source number density of 
UHECRs from the arrival distribution obtained by the AGASA. 

Figure 7 and figure 8 show simulated harmonic amplitudes. 
The number of simulated events is set to be 775 in the energy 
above 10 19 eV and their arrival direction is restricted in the 
range of -10° < 6 < 80° in order to compare our results with 
those of the AGASA. Note that 5 is the declination. The shaded 
regions represent 1 a total statistical error, which is caused by 
the two components of statistical error which occur from the 
finite number of simulated events and the random source se- 
lection from our IRAS sample. In order to see magnitudes of 
each error, we also draw errorbars, which represent the only 
statistical fluctuation due to the finite number of the simulated 
events. The event selection and the random source selection are 
performed 100 times and 40 times respectively. The regions 
below the solid lines are expected for the statistical fluctua- 
tion of isotropic source distribution with the chance probability 
larger than 10%. For all source number density, both first and 
second amplitudes show that our source models predict suffi- 
cient isotropy of UHECR arrival distribution obtained by the 
AGASA. 

Next, we investigate what number density of the sources re- 
produce the two point correlations obtained by the AGASA 
best. In order to evaluate it, we introduce xe m „ for a source 
distribution as 



X^max 



6=0 



{N(9)-N obs (0)Y 



(18) 



<d><e+A0 



where N(9) is the two point correlation function calculated 
from simulated arrival distribution within -10° < S < 80° and 
N o b s (0) is that obtained from the AGASA data at angle 9. a{9) 
is total statistical error of N(9) due to the finite number of sim- 
ulated events. The random event selection are performed 100 
times. This xo m „ represents goodness of fitting between the 
simulated two point correlation and the observed one. In this 
study, we take (9 max to be 10°. 
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log(E [eV]) 

FIG. 7. — The first amplitudes predicted by our source models, as a function of the cosmic ray energies. The shaded regions represent 1 a total statistical error 
due to the event selection and the source selection. The errorbars represent the statistical fluctuations only due to the event selection. The number of simulated 
events is set to be equal to that observed by the AGASA. The upper left is the first amplitude calculated by a source model whose number density of source Source 
is ~ 1 X 1CT 5 Mpc~ . The upper right, the lower left and the lower right are the first amplitudes for Source ~5x 1CT 6 , 2 X 10"*, 1 X 1CT 6 Mpc~ 3 respectively. 




Fig. 8.— 



The same as Fig. 7, but for amplitudes of the second harmonics. 
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10-» 10- 1 
Number density of sources (Mpc 3 ) 

FIG. 9. — xio as a function of number density of source « SO iirce)- The error- 
bars represent the statistical fluctuations due to the source selection from our 
galaxy sample. The left and right panel are calculated in the energy range of 
E > 4 X 10 19 eV and E > W [g eV respectively. 

Figure 9 shows xio as a function of number density of the 
sources (n SOU rce)- The errorbars represent the statistical fluctua- 
tions due to the random source selection from our galaxy sam- 
ple. The random source selections are performed 40 times. The 
two point correlation functions of the left and right panel are 
calculated in the energy range of E > 4 x 10 19 eV and E > 10 19 
eV, respectively. However, we cannot know the normalization 
of two point correlation function obtained from the AGASA 
for E > 10 19 eV. Thus, we fit the two point correlation func- 
tion obtained by the AGASA data to that calculated from our 
simulation at 8 m . dx . 

Source models with larger source number density have strong 
peak at a small angle scale on the two point correlation since 
there are some sources near to the earth. On the other hand, 
in the case of smaller source number density, a small number 
of sources near to the earth contribute the arrival distribution, 
especially in the highest energy case (left panel), since radial 
distances between any two sources from the earth are more dis- 
tant. Thus the peak of the two point correlation function also 
becomes more strong. Therefore xios should have a minimum 
as a function of source number density. As is seen from the 
figure, n source ~5x 10" 6 Mpc" 3 reproduces the two point corre- 
lation function obtained by the AGASA best. 

Therefore, Source ~5x 10~ 6 Mpc" 3 is the most appropri- 
ate number density of UHECR sources, since this source 
model also reproduces the harmonic amplitude obtained by the 
AGASA well. Note that number density of the sources have 
some uncertainty since the error bars in both panels are large. 

Figure 10 is the energy spectra at the earth predicted by 
our source models. These spectra are averaged ones among 
40 source distributions on each source number density. Solid 
line represents energy spectrum obtained for n source ~ 5 x 10~ 6 
Mpc" 3 . We also show the observed cosmic -ray spectrum by the 
AGASA (Hayashida et al. 2000). These simulated energy spec- 
tra have cutoffs around E ~ 10 19 6 " 8 eV except « SOU rce ~ 1 x 10" 5 
Mpc" 3 , which are the GZK cutoff. These spectra can repro- 
duce the AGASA data below 10 20 eV. Note that the spectrum 



with n source ~ 1 x 10" 5 Mpc" 3 has little spectral cutoff since 
there are many distributions which contain sources in the GZK 
sphere due to large source number density. The source model 
("source ~ 5 x 10" 6 Mpc" 3 ) also reproduces the observed energy 
spectrum only below 10 20 eV. We concluded in Yoshiguchi et 
al. (2003a) that a large fraction of cosmic rays above 10 20 eV 
observed by the AGASA might originate in the top-down sce- 
narios. Thus we consider UHECRs with only E < 10 20 eV in 
what follows. 



~1 xlO- 5 Mpc- 3 
~5 xl0-s Mpc- 3 
~2 XlO"" Mpc- 3 
~1 xlO- 6 Mpc- 3 
AGASA 




FIG. 10. — Energy spectra at the earth predicted by our source models. 
These spectra are averaged ones 40 times of the source selection on each 
source number density. Solid line represents energy spectrum obtained by a 
source model whose number density is ~ 5 X 10~ 6 Mpc -3 . We also show the 
observed cosmic-ray spectrum by the AGASA. 



5.2. Arrival Distribution of UHECRs above 10 19 eV 

In this subsection, we demonstrate a skymap of the arrival 
distribution of UHECRs in the case of n S0U rce ~5x 10~ 6 Mpc" 3 . 
We construct their arrival distribution using our method ex- 
plained in the section 4.3. 

Figure 1 1 shows one of results of the event generation above 
10 19 eV calculated from a specific source model with n source ~ 
5 x 10~ 6 Mpc" 3 . The points represent each event and the events 
are colored according to their energies. The number of events is 
5000, which is the expected number of events observed by the 
Pierre Auger observatory for a few years (Capelle et al. 1998). 
The skymap generated with both the EGMF and the GMF is in 
the lower right panel and that without any magnetic fields, that 
with only the EGMF and that with only the GMF is in the upper 
left, the upper right and the lower left panel respectively. 

This specific source model has three strong sources (see 
upper right). One is (l,b) ~ (199°, 34°), another is (l,b) ~ 
(287°, 19°) and the other is (l,b) ~ (25°, 11°). Each distance 
from us is about 77 Mpc, 65 Mpc and 70 Mpc respectively. In 
the absence of any magnetic fields (upper left panel), there are 
the strong clusterings of events at the directions of these three 
sources. When the effects of the EGMF are included (upper 
right), we find the diffusion of the clustered events. In the lower 
left panel, the clustered events are arranged in the order of their 
energies, reflecting the directions of the GMF. This was pointed 
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Without the EGMF and with the GMF 



With the EGMF and the GMF 



FIG. 1 1 . — Skymaps of arrival distribution of UHE protons with E > 10 19 eV at the earth, which is expected for the source model of uppermost panel in the 
galactic coordinate. We show only the sources within 200 Mpc from us for clarity as circles of radius inversely proportional to their distances. The events are shown 
by color points according to their energies. The event number is 5000, which is the expected number of events observed by the Pierre Auger observatory (Capelle 
et al. 1998) for a few years. The upper left, upper right, lower right and lower left panel is calculated without any magnetic fields, with only the EGMF, with only 
the GMF and with both magnetic fields respectively. 
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out by Alvarez-Muniz, Engel & Stanev (2002) and Yoshiguchi 
et al. (2003b). Note that we cannot find the clustered events at 
direction of one of the strong sources (l,b) = (25°, 11°). This 
is because UHE protons injected at this source cannot reach the 
earth due to the GME In the lower right panel, we also find the 
arrangements at the same points of the lower left panel. In ad- 
dition, the EGMF diffuses these clustered events as we see in 
the two upper panels. 

In order to see these features quantitatively, we compare the 
statistical quantities calculated with the EGMF to those calcu- 
lated without the EGMF in the presence of the GMF in the next 
subsection. 



5.3. Statistics on the UHECR Arrival Distribution 

In this subsection, we compare the statistical quantities on 
the arrival distribution calculated with the EGMF to those with- 
out the EGMF. We take n source - 5 x 10" 6 Mpc" 3 . 

Figure 12 and figure 13 show the two point correlation func- 
tions simulated by our source model in the energy range of 
E > 4x 10 19 eV and E > 10 19 eV respectively. In each figure, 
the left panel shows the two point correlation function calcu- 
lated with the EGMF and the right panel shows that without 
the EGMF. Note that the GMF is taken into account in the fig- 
ures. We calculate the two point correlation function for the 
simulated events within only -10° < S < 80° in order to com- 
pare our results with the AGASA data. The shaded regions 
represent 1 a total statistical error, which is caused by the finite 
number of simulated events and the random source selections 
from our IRAS sample. We also draw errorbars, which repre- 
sent the statistical fluctuations due to the finite number of the 
simulated events, which is set to be equal to that observed by 
the AGASA (49 events for E > 4 x 10 19 eV and 775 events for 
E > 10 eV). The event selection and the source selection are 
performed 100 times and 40 times respectively. The histograms 
represent the AGASA data (Takeda et al. 2001). For the data of 
E > 10 19 eV, we normalize the two point correlation function as 
the correlation function obtained by the AGASA fits the calcu- 
lated one at 30°, since we cannot know the normalization of the 
AGASA data with this energy. 

In both figure 12 and 13, it is visible that a peak at small 
angle is much stronger than that of the AGASA though the 
AGASA data were covered in the shaded regions, which are 
mainly caused by the source selection. In our previous work 
(Yoshiguchi et al. 2003b) in which the EGMF was neglected, 
we also faced this situation and pointed out possible explana- 
tions, one of which was effects of the EGMF. 

In figure 1 3, we find that a peak of two point correlation func- 
tion calculated with the EGMF at small angle scale is weaker 
than that without the EGMF. Thus the consistency with the 
AGASA data becomes better due to the EGMF. On the other 
hand, in figure 12, a obvious difference between the two panels 
on the peak at small angle is not found. This is because UHE- 
CRs above 4 x 10 19 eV are less deflected by the EGMF and 
hardly dispersed. 



1000 




6[deg] 



FIG. 12. — Two point correlation functions calculated from our source model 
with n s ~ 5 X lO^Mpc" 3 for E > 4 X 10 19 eV (49 events). The left panel 
shows the two point correlation calculated with the EGMF and the right one 
shows that calculated without the EGMF. The shaded regions represent 1 <r 
total error to consist of statistical error due to the finite number of the simu- 
lated events, which is set to be equal to that observed by the AGASA within 
-10° < S < 80°, and one due to the source selection. The errorbars represent 
the statistical fluctuation due to the finite number of the simulated events in 
order to see this error contribute to the total error. The histograms represent 
the AGASA data. 



1000 




10 20 10 20 



9[deg] 

FIG. 13. — Same as figure 12. But these two point correlation functions are 
calculated for E > 10 19 eV (775 events within -10° < 8 < 80° ). We fit the 
AGASA data to that calculated from our source model at S m ax = 30° since we 
cannot know the normalization of the AGASA data with this energy range. 

As mentioned above, the calculated two point correlation 
functions have large errors since some source distributions out 
of 40 contain very near source in our source model. Such source 
distributions do not reproduce the large-scale isotropy observed 
by the AGASA. We check that 20 source distributions out of 
40 predict the sufficient large-scale isotropy obtained by the 
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AGASA within 1 a statistical error due to the finite number 
of the simulated events. From these 20 source distributions, we 
calculate the two point correlation function in the presence of 
the EGMF again. The results are shown in figure 14. 



K > 4 xlO 19 eV 
4 9 events 



E > 10 19 eV 
775 events 



20 

Sfdee 



FIG. 14. — Two point correlation functions calculated with the EGMF from 
20 source distribution which satisfy the large-scale isotropy observed by the 
AGASA sufficiently, for E > 4 X 10 19 eV (49 events, left) and E > 10 19 eV 
(775 events, right). These functions are calculated within —10° < <5 < 80°. 
The shaded regions and the errorbars have the same mean as figure 12, 13. 
The histograms represent the AGASA data. The normalization of the AGASA 
data with E > 10 eV performs the same as figure 13. 

In figure 14, we find that the consistency with the AGASA 
data becomes better than the two point correlation functions in 
figure 12 and 13. We also see that the errors due to the source 
selection become small. However, we should note that the peak 
at small angle scale is still relatively strong, compared with the 
AGASA though our previous result (Yoshiguchi et al. 2003b) 
is improved by the effect of the EGMF. We assume effects of 
the random component of the GMF, which is neglected in this 
work, as one of possible explanations for this fact. This issue is 
left for future investigations. 

We also investigate the harmonic amplitudes in the same 
way. But there is little difference dependent on the EGMF. 

6. SUMMARY AND DISCUSSION 

We presented numerical simulations on propagation of UHE- 
CRs above 10 19 eV in a structured EGMF and GMF. We used 
the IRAS PSCz catalogue in order to construct a structured 
EGMF model which reflects the local structures actually ob- 
served, and source models of UHECRs. We calculated an 
inverse process of their propagation taking the energy loss 
processes into account in the EGMF. We injected UHECRs 
from the earth isotropically whose charges are taken as -1 and 
recorded these trajectories. They could be regarded as trajec- 



tories of UHE protons from the extragalactic space. We then 
select some of their trajectories according to given source dis- 
tributions. The simulated arrival distribution was able to be ob- 
tained by mapping the velocity directions of the selected tra- 
jectories at the earth. The use of this method enabled us to 
calculate only trajectories of UHECRs reaching the earth and 
saved the CPU time effectively. The validity of this method 
was supported by the Liouville's theorem. 

We calculated the harmonic amplitudes and the two point 
correlation functions of arrival distribution of UHECRs above 
10 19 eV, using our source models and examined what number 
density of the sources reproduces the large-scale isotropy and 
the small-scale anisotropy obtained by the AGASA best. As a 
result, we found that ~5x 10~ 6 Mpc -3 was the most appropri- 
ate number density of source of UHECRs. Number density of 
source is a constraint on source candidate of UHECRs. 

We also demonstrated skymaps of the arrival distribution of 
UHECRs above 10 19 eV, using our source model for « SO urce ~ 
5 x 10~ 6 Mpc" 3 with the event number expected by future ex- 
periments and examined how the EGMF affects the arrival dis- 
tribution of UHECRs. The main result was diffusion of clus- 
tering events which is obtained by calculations neglecting the 
EGMF. 

In order to see the effect of the EGMF quantitatively, we 
compared the two point correlation function calculated with our 
structured EGMF model to that without the EGMF. We found 
that the EGMF weakened the small-scale anisotropy and im- 
proved a prediction in Yoshiguchi et al. (2003b), which had 
been calculated with only the GMF. 

However, we found the calculated two point correlation func- 
tions had large errorbars since source distributions, which con- 
tained sources very near to us, existed. Such source distribu- 
tions do not reproduce the large-scale isotropy observed by the 
AGASA. Thus we calculated the two point correlation func- 
tions from source distributions which predicted the large-scale 
isotropy obtained by the AGASA within 1 a statistical error. 
These simulated two point correlation function reproduced that 
of the AGASA better. It is possible that these functions at small 
angle scale can be closer to the observational data due to the 
random component of the GMF. This issue is left for future 
studies. 

New large aperture detectors are under development, such as 
the Pierre Auger observatory (Capelle et al. 1998), the EUSO 
(Benson & Linsley 1982) and the Telescope Array. These 
projects are expected to increase observed events of UHECRs 
per year drastically. We can obtain more strong constraints on 
our source model, other than number density of sources, using 
other statistical quantities when the detailed data of large events 
of UHECRs are published. This is one of plans of future stud- 
ies. 
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